Speed of evolution in large asexual populations with diminishing returns 
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The adaptive evolution ol large asexual populations is generally characterized by competition 
between clones carrying different beneficial mutations. This interference phenomenon slows down the 
adaptation speed and makes the theoretical description of the dynamics more complex with respect 
to the successional occurrence and fixation of beneficial mutations typical of small populations. A 
simplified modeling framework considering multiple beneficial mutations with equal and constant 
fitness advantage captures some of the essential features of the actual complex dynamics, and some 
key predictions from this model are verified in laboratory evolution experiments. However, in these 
experiments the relative advantage of a beneficial mutation is generally dependent on the genetic 
background. In particular, the general pattern is that, as mutations in different loci accumulate, 
the relative advantage of new mutations decreases, trend often referred to as "diminishing return" 
epistasis. In this paper, we propose a phenomenological model that generalizes the fixed-advantage 
framework to include in a simple way this feature. To evaluate the quantitative consequences of 
diminishing returns on the evolutionary dynamics, we approach the model analytically as well as 
with direct simulation. Finally, we show how the model parameters can be matched with data from 
evolutionary experiments in order to infer the mean effect of epistasis and derive order-of-magnitude 
estimates of the rate of beneficial mutations. Applying this procedure to two experimental data sets 
gives values of the beneficial mutation rate within the range of previous measurements. 



I. INTRODUCTION 

Although the links between the statistical theory of 
evolution and statistical physics can hardly be consid- 
ered novel ^ , in the recent years they are attracting the 
attention of researchers in a renewed way. What changed 
is the increased availability of genomic and experimental 
data, in amounts and levels of precision that (despite 
there is still a considerable room for improvement), were 
difficult to envisage in the past century. These innova- 
tions are consolidating the biology and genomics of evo- 
lution into a more stable field of application of modeling 
ideas from statistical physics. 

In particular, contemporary technology allows the re- 
alization of controlled laboratory evolution experiments, 
which can guide theoretical investigations and in princi- 
ple makes the validation and falsification of phenomeno- 
logical theories feasible [2- These experiments are often 
performed with large populations of microorganisms, and 
allow to explore adaptation under well-defined sources 
of natural selection. Moreover, phenotypic characteri- 
zation and high-throughput sequencing give a quantita- 
tive insight into the genomic adaptation of these micro- 
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bial populations [31 13], with notable consequences in a 
wide range of bio-technological and ecological contexts. 
For these large asexual (or rarely mating) populations, a 
high number of beneficial mutations emerge in different 
clones, and cannot be mixed because of slow or absent 
recombination. These beneficial mutations appearing in 
parallel coexist and compete to drive adaptation. This 
phenomenon of concurrent beneficial mutations (some- 
times generically termed "clonal interference" ) , is related 
to the Fisher-Muller hypothesis (or Hill- Robertson effect) 
for the advantage of recombination . In general, ben- 
eficial mutations also arise with a distribution of fitness 
advantage, which is generally believed to be exponen- 
tial [S]. 

Recent models have generally dealt with the compe- 
tition between mutations of different strengths and the 
competition between mutations that arise on different fit- 
ness backgrounds separately. The first effect, the role of 
a distribution of fitness effects, is analyzed by so-called 
"clonal interference" models [37] [TJ [71 [H] . In these mod- 
els any individual is either the wild type or a mutant 
derived directly from the wild type. Thus, multiple mu- 
tations arising from the extant clones are neglected. Con- 
versely, models that explicitly deal with multiple muta- 
tions typically assume that all mutations have the same 
effect [1 , 9 1 Ij . The latter kind of model has the advan- 
tage of being simpler to treat and accessible analytically, 
and is characterized by a Gaussian-like traveling wave for 
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the histogram of log-fitness throughout the population. 
In absence of epistasis, this wave moves towards higher 
log-fitness with a constant speed and shape. New mu- 
tations are fixed in the population if they occur in the 
high-fitness edge (or "nose") of the distribution. Conse- 
quently, a quantitative law relates the width of the log- 
fitness histogram and the adaptation speed. Recent work 
incorporating both effects [T^] has shown that if the dis- 
tribution of fitness advantages given by beneficial muta- 
tions is sufficiently peaked around a single characteristic 
value, the evolutionary dynamics can be described by an 
effective theory with a single typical selection coefficient 
and a rescaled beneficial mutation rate. 

While the models described so far have been matched 
successfully with the diversity and adaptation speed of 
short-time laboratory evolution experiments [I3J, there 
appears to be one important discrepancy between the 
models and the behavior of bacteria evolved in the labo- 
ratory for longer times (roughly, > 1000 generations). 
Furthermore, two recent experimental studies [141 115] 
have shown a general pattern in the advantage of com- 
bined beneficial mutations occurring in different genes. 
This combined advantage is lower than the sum of that 
of individual mutations. In other words, when muta- 
tions of loci in different genes accumulate, the effective 
advantage of each of them is lower. This was shown by 
combinatorial genetics techniques, by constructing all the 
possible configurations of a small set of mutations, and 
evaluating their advantage through competition experi- 
ments. This trend, referred to as "diminishing returns" 
epistasis, had been previously suggested theoretically on 
the basis of the general pattern of adaptation observed in 
long-term microbial experiments jl6j . using a modeling 
framework that neglected concurrent or multiple muta- 
tions. Another study predicts the same principle on the 
basis of a simple fitness landscape model combined with 
the distribution of single mutation effects measured ex- 
perimentally [17j . The actual pattern in the fitness as- 
sociated to the same mutation in different backgrounds 
observed by the two studies is complex, as, on top of 
the diminishing return effect, the advantage appears to 
depend on the mutation identity. Based on their data, 
Chou and coworkers [H] defined an expression for the 
fitness made of two additive components. They consider 
the growth rate as the sum of two "constituent pheno- 
types" : metabolic rate and protein expression burden. 
When mutations are combined, the metabolic and ex- 
pression contributions to the fitness combine multiplica- 
tively in an independent way. They show that this model 
generates an accurate prediction for the fitness of all their 
multi- allele strains. Even more recent systematic experi- 
ments [4 are unveiling a complex scenario where different 
mechanisms coexist for the interactions of mutations be- 
tween and within functional "blocks" , which can span 
multiple genes along the genome. One interesting mod- 
eling approach that was developed recently [18] incorpo- 
rates genetically linked multiple mutations explicitly, and 
describes interference interactions between multiple ben- 



eficial and deleterious mutations. This description has 
enough degrees of freedom to possibly account for many 
of the phenomena observed experimentally. However, the 
full experimental complexity is difficult to incorporate in 
a treatable model, and experimental data on linked mu- 
tations and interference between them arc not easy to 
be obtained. A simplified descriptions in the spirit of 
the multiple mutations model can be useful in order to 
model evolving populations using a minimal quantity of 
information on mutations and fitness advantage. 

Here, we take a simplified approach to investigate the 
diversity and speed of adaptation in presence of dimin- 
ishing returns. We define a framework that can account 
for multiple mutations, and incorporates the effect of di- 
minishing return epistasis in a simplified way. Namely, 
the fitness of a mutation depends only on its order of 
appearance in a clone, and decreases with it. This gen- 
eralizes the standard multiple mutation model, which we 
recover in the case in which the fitness decrease with mu- 
tation index is zero. We preserve the model assumption 
that evolution is driven by beneficial mutations which ap- 
pear with constant rate (see the Discussion for an eval- 
uation of these assumptions in light of the results) . We 
study the mean-field behavior of this model using stan- 
dard techniques [T], and derive a relation that is analo- 
gous to Fisher's theorem. Moreover, we discuss how the 
model assumptions cease to be valid when the rate of 
fixed beneficial mutations drops to the level where dele- 
terious mutations should be accounted for. We also ex- 
plore the behavior at finite population sizes, and obtain 
simple analytical quantitative estimates of the speed of 
adaptation that generalize that of the standard multi- 
ple mutation model. Finally, we show how the model 
parameters can be matched with data from two differ- 
ent long-time laboratory evolution experiments and how, 
assuming the model, one can derive order-of-magnitude 
estimates of the beneficial mutation rate and the mean 
effect of the diminishing returns. 



II. BASIC FEATURES OF THE MODEL 

A. Model definition 

We model an asexual population of N haploid indi- 
viduals, or sequences, in which each individual of type 
i produces a random number of offspring with average 
equal to its fitness Wi. Inheritance is introduced by as- 
signing the fitness of their parents to the offspring. The 
evolutionary dynamics of the population is based on the 
well-known Wright-Fisher model [111 [20] • We used the 
following reproduction scheme |21j . Each individual at 
generation 1 is chosen as the offspring of an individual 
in class i present at generation t with probability Xi/-^i 
where Xi — Wi/{w) is the relative fitness of class i. This 
prescription for the evolutionary dynamics assumes non- 
overlapping generations since at each generation there is 
a complete replacement of parents with progeny. The 



3 



L(k) 



..llllll 




Number of mutations k 



Fitness Advantage Sag(k) 



B 



Ml. 



■ n 



ft) 



If 

lit 



Establishment 



'^ma/f+^y size -J/ZiS 



Number of mutations k 



Beneficial 
mutation rateL/^ 
Fitness w,=e-°" 


|Frequency f(k,t)\ 
1 


Surviving offspring 




probabilities p(k,t+l) 




1 



f(k,t+l) from 
multinomial distribution 



FIG. 1: (Color online) Basic features of the model. 

(A) Because of competition between beneficial mutations, the 
population is divided into sub-populations with different fre- 
quencies (left panel), defined by the number of mutations k. 
Lfc is the difference, in number of mutations, between the 
maximum number of mutations found in a clone femax and the 
mean. This induces a distribution for the log-fitness sog{k) 
(right panel). Both distributions travel in time, driven by es- 
tablished beneficial mutations, with instantaneous velocities 
Vk and Vs. (B) The dynamics is driven by beneficial mu- 
tations and the rate of establishment of new classes. The 
dashed rectangles represent the class frequency histogram at 
a subsequent time. (C) Sketch of the algorithm used in the 
simulations (from ref. [I], see text). The instantaneous fre- 
quencies of mutation classes define, through Eq. [s] (which in- 
corporates selection and mutation) the probability p{k, t + 1) 
that an individual with k mutations is found at the subse- 
quent generation. The frequencies at the subsequent genera- 
tions are then sampled from a multinomial distribution with 
parameters p{k,t + 1). This procedure allows to access large 
population sizes. 



population size N is kept constant, and there is no re- 
combination move. The presence of fitness differences 
in the population leads to natural selection since classes 
of individuals with higher fitness will generate increas- 
ingly larger fractions of the population as the dynamics 
proceeds, while classes with low fitness will progressively 
disappear. 

Each offspring has a constant probability per genera- 
tion Ub (the "beneficial mutation rate") of acquiring a 
beneficial mutation. Beneficial mutations increase the 
parental fitness Wi to a higher value w[ following the re- 
lation w'^ — Wi{l + s) ~ WiC^, where the positive param- 
eter s is the "selection coefficient" whose typical values 
in laboratory evolution experiments with bacteria are in 
the range s ~ 0.001-0.005 [H US]. We assume that each 
new beneficial mutations hits a new site on the genome. 



While the fitness of new mutations is a complex issue [23] , 
in presence of abundant beneficial mutations, deleterious 
mutations (negative effect on fitness, s < 0) do not typ- 
ically contribute to the adaptation of large populations 
and are customarily neglected [1] [TUl E] . 

The conditions for the emergence of the interference 
phenomenon between mutations for a large population 
can be understood with simple scaling arguments as a 
competition between processes occurring at different time 
scales p^. When a new mutant with advantage s arises, 
the probability that its lineage grows sufficiently in size to 
overcome genetic drift (stochastic fiuctuations in the re- 
productive process) and start to deterministically expand 
in the population (i.e. it is "established") is tt{s) ~ cs, 
where c is a constant factor that depends in the specific 
model used [25]. Fo r the algorithm used here, n{s) ~ 2s 
(see Appendix, Fig. A. 2) [3S]. Furthermore, since the ex- 
tinction probability is 1 — 7r(s), one can also estimate the 
population size conditioned to certainty of survival of the 
lineage as I/tt ~ l/2s. In absence of additional beneficial 
mutations, once a lineage is established (i.e. has survived 
genetic drift and has roughly size l/2s) it will take over 
the population (go to "fixation") logistically at initial 
rate s. Consequently, the scale of the fixation time can 
be estimated by imposing that ^ exp(srfix) — N, giving 



a characteristic time 



ln{2Ns) 



[TU] . On the other 



hand, the time scale for appearance and establishment of 
a new beneficial mutation is Test — t^ttt^^ztttt . Therefore, 
when NUh ^-^ ^ 



NUb-!T{s) 

2in(2Ns) ("^fix < Tost) a beneficial mutation 
can fix before any new mutation can establish, making 
the evolutionary dynamics driven by successive sweeps of 
new lineages arising on an essentially clonal population. 
This regime is called "selective sweeps" or "periodic selec- 
tion" . Instead, a sufficiently large population with high 
beneficial mutation rate evolves in the opposite regime, 
in which multiple mutations can establish before fixation 
of any of them and interfere with each other. This is the 
regime considered here, which is believed to be relevant 
for laboratory evolution experiments with microorgan- 
isms. 

As mentioned in the introduction, the standard ap- 
proach of multiple mutation models is to consider con- 
stant advantage s = sq for all mutations. This en- 
tails that, in absence of epistasis, the fitness s^k can 
be assigned to individuals (sequences) with k mutations, 
which (roughly) propagate with rate l + soik— (k)) |11| . 
(k) is the average number of accumulated mutations in 
a single realization of the process, i.e. the average over 
the distribution of k shown in Fig. [l] For a range of pa- 
rameter values, this constant advantage framework can 
be seen as an effective theory for a model with a distri- 
bution of fitness for beneficial mutations [T^] . 

Following the spirit of this modeling framework, we 
build a minimal phenomenological model including di- 
minishing return epistasis in presence of competition be- 
tween beneficial mutations. The model assumes that suc- 
cessive mutations do not lead to the same fitness gain. 
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but the fitness gain is dependent on the number of mu- 
tations already occurred in an individual |16j . A sim- 
ple way to implement this feature is to extend the stan- 
dard model for multiple mutations and consider selection 
coefficients dependent on the number of mutations, i.e. 
s — sog'{k), where g'{k) is a decreasing function of k. 
As in the standard multiple- mutation model [101 E] > the 
population can be divided into classes of individuals with 
the same number of mutations that are in one-to-one cor- 
respondence to fitness advantage classes, as represented 
in Fig.[T]A_. An individual with k beneficial mutations has 
fitness 



beneficial mutations (see Appendix [A] for the derivation) 
has the form 



w, = e''>' . (2) 

Appendix |A] describes two additional examples of the 
diminishing returns function g'(fc) that we have consid- 
ered, leading to logarithmic (i.e. the case a = 0) or 
geometric increase of the log-fitness. The comparison be- 
tween the three variants will be useful in the second part 
of this work, when the model parameters are matched 
with experimental data. 



Since the rules defining the dynamics contain the relative 
fitness Xki the model is unaffected by multiplication of 
all the Wk by a common factor. Therefore, the fitness 
value wq — I can be arbitrarily assigned to the genotype 
with no mutation ("ancestral" or "wild- type" ) . 

We considered typical values of the parameter sq be- 
tween 10""^ and 10~^ [ll[23]. For the parameter Ub, we 
explored the range 10~^° — 10~^ [HI [53], and we con- 
sidered population sizes between 10^ and 10^° [1]. These 
ranges, which can be considered plausible with respect to 
what is known empirically, impose a hierarchy of scales in 
the effective parameters. For example, typically Uh <^ Sq 
and sqN is large. As for UbN, we are interested in 
exploring sufficiently large values to ensure the regime 
of concurrent mutations. However, since the advantage 
given by a beneficial mutation SQg'{k) decreases during 
the time-evolution of the model, this hierarchy can in 
part change over time. When the advantage becomes 
too small, i.e. for large times, some of the assumptions 
of the model become unrealistic. Namely, for sufficiently 
large times, deleterious mutations cannot be neglected, 
and the establishment size l/7r(s) can eventually become 
comparable to TV. The inequality NUt ^ 2in(2Ns) ' "^^lich 
ensures concurrent mutations, can also break down for 
very large times, when s decreases. This is not a prob- 
lem - as any experiment spans a finite time. However, in 
contrast to the standard multiple mutation model, which 
describes a stationary state where the regime is entirely 
defined by the parameters, in this model the time-range 
of validity of the assumptions has to be kept in mind. 
We will comment on these issues in the following sec- 
tions, when reviewing the results. 

In order to fully define the model, one has to choose a 
specific functional form for the function g'{k), describing 
the strength of the negative epistasis between mutations. 
In general, every function g'(fc) leading to a sum g{k) — 

Sfe'=o ^09' {k') that is sub- linear in k defines a model with 
diminishing return. A simple example is given by the 
choice of a fitness gain that depends on the number of the 
extant mutations as a power law. In this case, g'{k) — 
ak"~^ with a < 1, where the epistasis grows in strength 
with decreasing a from the non-epistatic case a = 1. In 
the case < a < 1, the fitness of an individual with k 



B. Background on the multiple mutation model 
with no epistasis 

For a = 1, the power law return model reduces to the 
particular case of absence of epistasis (i.e. g'{k) — 1, 
hence = e^^°). This case coincides with the stan- 
dard multiple mutation model for the evolutionary dy- 
namics of a large asexual population [TUl [2^ , in which 
every mutation carries a constant fitness advantage sq- 
The phenomenology of the multiple mutation model with 
constant advantage has been explored in detail theoreti- 
cally [ini [U and some of the resulting predictions have 
been experimentally tested |13| . 

As mentioned in the Introduction, one of the main re- 
sults of these studies is that a steady state exists, where 
the population is organized in fitness classes (and the cor- 
responding mutation classes, as illustrated in Fig. [l]A_), 
forming an approximately Gaussian traveling wave en- 
compassing a constant number 2L of mutation-classes 
and propagating to higher average fitness with constant 
velocity Vs- Technically, this wave has all the character- 
istics of a soliton [T] [2S] . Importantly, the wave is driven 
by the mutations that are fixed at the edge of the distri- 
bution of the advantage. Fig. [Tj3 depicts schematically 
how the fitness and mutation class waves propagate. Mu- 
tation class k is fed by beneficial mutations from the pre- 
vious class, k — \. The fittest class kmax produces a new 
class of mutants that will be subjected to genetic drift 
until its size reaches the establishment size that depends 
on its relative advantage. Assuming new class kmax + 1 
senses a constant background from the class with mean 
fitness, the advantage can be estimated as SqL, where L 
indicates its distance from {k) (see Fig. [T]). Under this 
argument, the establishment size is l/2As ^ 1/2sqL. 
During the time in which a new fittest class establishes, 
the peak of the distribution will advance and the less- 
fit clones die out, thus keeping constant the distribution 
width (i.e. constant L). 

The speed of adaptation can be estimated by a scaling 
argument [TJ [TOl [H] imposing that the time r during 
which a new class establishes (and the class histogram 
moves by one bin) has to correspond to the larger-scale 
movement of the histogram by deterministic growth. 
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The first part of tlie argument estimates t from the 
condition 



1 



dt e 

2soL 



so{L-l)t 



tt{soL) 



(3) 



which imposes that order one mutations are produced 
(and establish with probability tt) from the foremost bin 
while it grows exponentially according to its advantage. 

The foremost bin is treated as it were growing in a 
constant background with mean fitness (as above) and 
1 /2sqL is taken as the size of the foremost bin at its birth. 
Integration of the above equation yields the expression 



1 SqL 



which is valid for sufhciently large L and advantage of 
the fittest class, since the contribution of the integration 
boundary t = has been neglected. 

The second part of the argument can be formulated as 
a normalization condition of the fitness class histogram, 
joint with the condition that each bin grows exponen- 
tially according to its advantage, 

^ ^ ^ f 1 I e'*o(i-l)r , so((L-l) + (L-2))r , \ 

2 2soL V 7 

note that this condition neglects the contribution of 
beneficial mutations from the neighboring bins. Approx- 
imating the above formula as an integral computed with 
saddle-point (i.e. assuming that the largest term in the 
sum, the class with mean fitness, dominates [101 HI]) 
yields the expression 



log{soLN) 



2 



(6) 



Equating the two expressions Eq. Q (or some variants 
that can be derived more rigorously) and ([6) for L or t 
generates implicit estimates for the speed of the mutation 
classes histogram, Vk — 1/t, which in general work very 
well, despite of the approximations taken. In particular, 
neglecting the addends logL in Eq. Q and (|6), i.e. up 
to logarithmic corrections in L, one obtains the closed 
formula fTUl H] 



]_ ^ (log so/Ub)^ 
Vk 2so log(iVso) 



(7) 



Finally, Vg can be estimated by sq/t. 

In the general case of diminishing return epistasis, 
the variables of interest become a function of the ge- 
netic background. Specifically, the speed of adaptation 
in log-fitness space Vs and the speed Vk of the mutation 
classes, are two non-trivially distinct quantities, and the 
same holds for the width of the mutation and the advan- 
tage histograms {Lk and Lg) relative to the population 
(Fig.[l]). 



C. Simulation algorithm and effective parameters 

The model was simulated with the algorithm of Park 
and Krug, as described in refs [U [5T]. The simulation 
scheme is sketched in Fig. [Tp. In a typical initial con- 
figuration, all clones have k = mutations. At each 
subsequent time step the progeny of the individuals of 
all fitness classes are sampled from a multinomial distri- 
bution of parameters {p{kjt + l)}{fee[fe,„i„,fe„ax]}- The pa- 
rameters p{k, t + 1) take into account the relative fitness 
Xk computed at time t and the contribution of beneficial 
mutations. Specifically 



(4) pik,t + l) = {1-Ub)fik,t)xkit) + Ubf{k~l,t)xk- 



(8) 



(see also Eq. ([9|)). Together, these definitions are equiv- 
alent to a Wright-Fisher model with separate selection 
and mutation steps. In particular, the microscopic step 
described above reduces to the standard Wright-Fisher 
model for genetic drift when Ub — 0. The multinomial 
random numbers are generated by iteratively drawing bi- 
nomial random numbers with parameters q{k, t + 1) — 
p{k, t+l)/ J2k Pi^^ t + 1) starting from fcmax- 

The model defined above is invariant by suitable rescal- 
ing of time, provided the other model parameters are also 
rescaled correctly. This feature is useful in comparison 
with experiments (see following) in order to understand 
the correspondence between time steps in the model and 
a generation in the experiment. It can also be useful in 
increasing the efficiency of simulations. 

Let us suppose one desires to map a reference empiri- 
cal time into the time steps of the model. This requires 
the mapping tmodei = i^temp, in the expression describ- 
ing the growth of a sub-colony from a clone in terms 
of the experimental generation time. The coefficient r 
is the number of empirical generations corresponding to 
a time step, and represents the simulation time scale. 
A simple choice of time scale is the case r — 1, imply- 
ing a one-to-one correspondence between empirical gen- 
erations and time steps in the model. Since an estab- 
lished clone with advantage s grows as e''*, the advan- 
tage function is proportional to the time scale as follows, 
Smodci = f Scmp- This means that the value of the fit- 
ness in the model depends exponentially on r. Similarly, 
since the beneficial mutation rate is defined as the num- 
ber of expected mutations per genome per generation, 
the map between time steps and empirical generations 
implies C4,cmp = model/''- Finally, the correct rescaled 
population size is A^modci = A'emp/?'- In practice N can 
be rather large in a typical experiment (e.g. ~ 10^), so 
that the rescaling of population size does not affect much 
the dynamics provided r is not too large. The rescaling 
of the parameters described above can easily be ratio- 
nalized keeping in mind that the basic time scales of the 
model are set by the products sN and NUb- In order to 
verify that the invariance discussed above is valid for our 
model, where the advantage s varies with the number of 
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mutations as sog(fc), we ran some simulations choosing 
the model parameters {N, sq, Ub) using diff erent maps 
between the time units (see Appendix, Fig. A.f ). The 
results indicate that for any practical purpose, the in- 
variance is effective. 



III. RESULTS 



Adaptation slows down until it reaches an 
effective arrest. 
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FIG. 2: (Color online) The mean number of mutations 
and the mean advantage grow sub-linearly with time. 

The top panels show the increase in time of {sog{k)) (A) and 
of (k) (B) obtained by direct simulation of the diminishing 
return model. The plot in panel B is in log-log scale, and the 
data are compared with a reference line (dashed blue line, 
with slope 5 • 10~^) to highlight the sublinear growth of (k). 
The continuous red line shows the the asymptotic long-time 
linear behavior with slope corresponding to Ut- (C) Long- 
time behavior of the mean speed of fixed mutations vi^ (green 
symbols), averaged over difi^erent realizations. For long times, 
this quantity decreases (as a power law) towards the limit 
value Vk = Ub, where the assumptions of the model break 
down and deleterious mutations need to be accounted for [27j . 
This limit also corresponds to the limit value of obtained 
by a mean-field estimate (see text). Simulations are carried 
out using the parameters N — 5 ■ lO", so — 0.5, a — 0.02, 
Ub = 1- 10"^. Avera ges are computed over 100 realizations 
(these averages are implied in the notations for the y-axis 
labels) . 



the advantage. The average number of mutations (fc) also 
grows sub-linearly with time, for intermediate to long 
times (Fig. [2j3). This trend is independent from a (or 
from the specific model of the of decreasing g'{k) and is 
due to the fact that decreasing advantage and the con- 
sequent rise of the establishment threshold for clones to- 
gether slow down adaptation. The time derivatives of 
(sogik)) and (fc) estimate the adaptation speed Vs and 
the mutation-accumulation speed Vk of a typical realiza- 
tion. Figure [2p shows an average of Vk over one hundred 
realizations, plotted as a function of time. The simula- 
tions indicate that Vk relaxes to a plateau which is close 
to the beneficial mutation rate Ub- Equivalently, for long 
times, the mean number of fixed mutations shows a lin- 
ear behavior in time with a rate close to Uh (red line in 
Fig. [2J3). In the same long-time limit, the advantage of a 
mutation s — SQg'{k) drops asymptotically to zero. 

The long-time regime where (k) follows a linear trend 
is outside of the limit of validity of the model, and has 
to be regarded as unphysical, since when Vk = Ub, dele- 
terious mutations cannot be neglected [17]. Thus, for 
any finite N the asymptotic trend of (k) has to be inter- 
preted as an effective signature of an evolutionary arrest 
for both Vk and Vs , where beneficial mutations should be 
in equilibrium with deleterious one, which could possibly 
be captured by a variant of the model including dele- 
terious mutations [27]. Note also that realistically the 
beneficial mutation rate itself could decrease in the later 
stages of evolution [TJ [55] . 

Despite the fact that the assumptions of the model 
break down asymptotically, the long-time phenomenol- 
ogy appears to be theoretically consistent. Indeed, the 
long-time regime can be easily rationalized as a transi- 
tion to an "effectively neutral" evolution regime. In fact, 
the selection coefficient sog'{k) becomes irrelevant at long 
times, so that the fixation dynamics is driven solely by 
genetic drift. Given that deleterious mutations are not 
accounted for, and the benefit carried by the allowed mu- 
tations becomes negligible, in this regime Ub effectively 
becomes a neutral mutation rate. The probability of fix- 
ation of an essentially neutral mutation is ~ while 
the rate of appearance of new mutations is UbN. There- 
fore, the pace at which new mutations are accumulated is 
approximately Vk Ub- As we will see, this asymptotic 
velocity Ub is also recovered from a mean-field estimate 
of Vk, valid in the infinite population limit. 



B. A mean- field analysis gives a relation between 
the variance of the advantage distribution and the 
adaptation speed. 



Direct simulation of the model (Fig. [2]A) , shows that 
the mean advantage (sog(fc)) grows sub-linearly with 
time, as expected from the diminishing returns pattern of 



In the limit of infinite population, N 00 the dy- 
namics of the model can be described using the following 
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FIG. 3: (Color online) The mean field approximation 
captures a relation between Vk and the width of 
the fitness class distribution, which is valid at in- 
termediate times for moderate A'^ and until longer 
times for large population sizes. (A) Simulated vari- 
ance of the mutation classes distribution, shown as a func- 
tion of the expected variance from the mean-field estimate 
(VarA/F = {vk - Ub){ky'°'{soa)-\ see Eq ([ll|). To avoid 



ambiguities, averages over realizations are indicated by a suf- 
fix R. The continuous red line represents the theoretical pre- 
diction Vari\/F = Varfc. The error bars (standard deviations 
over realizations, anCVaxk)) become larger with (Varfc)_R. (B) 
While for increasing times (Tij(Varfc) diverges, the relative 
variability crj{(Varfc)/(Vari,)_R over realizations decreases with 
increasing population size A'^, for any fixed time. This sug- 
gests that the infinite A'' mean- field estimate is well-defined. 
Simulations are carried out using the parameters sq = 0.5, 
a = 0.02, Ub = 1-10"^. Population size in panel A is = lO'^. 



mean- field equation [T], 



Ubf{k-l,t-l) 



w(k) 
w)(t-i) 
w{k - 1) 

H(t-i) 



(9) 



where f{k,t) is the frequency of individuals with k ben- 
eficial mutations at generation t, w{k) — e""^'-'^-', and 
{'w)t = Y^kWkf{k,t) the mean fitness. 

Multiplying Eq. ([9]) by A; and summing over k gives 
the following expression for the dynamics of the mean 
number of mutations (fc)t = kf{k,t), 



{k)it+l) 



{w){t) 



(10) 



This expression can be further simplified assuming that 
the frequency distribution is narrowly peaked around the 
mean (which travels in time) and that it can be expressed 
as f{k,t) ss 6{k; {k){t)). In this case it can be easily 
verified that Vk = Ub- 

A different, more instructive, relation, which keeps into 
account the width of /(fc, t) can be obtained starting from 
and expanding w{k) under the assumption that 
{k — (fc)) <C (fc), for every index k of non-empty 
classes. This assumption is verified by simulations (see 



Eq. 



10 



Fig. A. 3 in Appendix) and by further considerations on 
the finite- width of the distribution given in the follow- 
ing sections. 



Computing averages in Eq. (10 1 to first order in Dk, 
and noticing that {Dk) = and (S^.) = Var^, it is possi- 
ble to obtain an expression for the speed of accumulated 
mutations as a function of the variance of the fitness class 
distribution. Estimating Vk as « {k){t -I- 1) — {k){t) 
gives, for the case g{k) — k" 



d{k) 
dt 



saa{k)'^-\t)\a.Vk{t) + Ub 



(11) 



According to this equation, Vk is driven by two terms, the 
increase of mutations due to the beneficial mutation rate 
Ub and the selection of individuals with larger fitness. 
This result is somewhat reminiscent of Fisher's funda- 
mental theorem and the Guess relation, which relate the 
speed of adaptation to the variance of the fitness. In this 
case, the speed of accumulation of successive mutations is 
related to the width of the mutation class histogram, but 
rescaled by the factor (fc"^^)(t), which decreases with 
time. For the non-epistatic case (a = 1) one recovers the 
more usual linear proportionality, since the advantage is 
linear in the mutation class index [1]. 

The mean-field limit of the model with a = 1 has been 
previously addressed by Park et al [T] with a moment 
generating function approach. In particular, they esti- 
mated the distribution variance as Vari, ~ i:i!Zt which 

So ' 

substituted in Eq. [IT] gives precisely their expression for 
the speed f ^ — 1 (in the limit of small Ub) suggesting that 
our result is a consistent generalization to the epistatic 
case. 

Therefore, for a — 1 the speed of adaptation does not 
depend on the mutation rate for infinite populations. On 
the other hand, for diminishing returns, the increase in 
the width of the distribution of fc, does not compensate 
for the term (which tends to 0), and, for long 

times, Eq. ( |11[ ) predicts the limit velocity Vk = Ub, as 
observed in simulations (Fig. [2]C). Since the increase in 
the width of the distribution depends on the size of the 
population and on the advantage sag{k), the unphysical 
limit velocity Vk ~ Ub is approached with different laws 
depending on the specification of g{k). 

It is necessary to discuss under which conditions this 
mean-field limit can be considered a valid estimate of the 
typical behavior of a realization. At finite N, simulated 
data are in good accordance for intermediate time with 
the mean field estimate for the variance Varj, of the mu- 
tation class distribution (Fig. [3|. Additionally, Varfe/(/c) 
decreases quickly with time (see Fig. A. 3 in Appendix), 
justifying the assumption of small Dk/{k) (since it is ex- 
pected that \Dk\ ^ (Varfe)i/2 for every k). However, 
the variability of Var^ over different realizations increases 
quickly. In order to verify whether these fiuctuations are 
well-behaved, we have evaluated 7 = an (Var^) /(Varfc)^, 
where Var^. indicates the variance of the mutation class 
distribution in a single realization (roughly analogous to 
L^), and the suffix R indicates averages over realizations. 
Specifically, {x)ii indicates the average of the quantity x 
over different realizations, while (7ji{x) is its standard 
deviation. Therefore, 7, plotted in Fig. |3l3 as a function 
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of N, represents the relative variability over the realiza- 
tions of the variance Var^ of the distribution. For any 
fixed time, this quantity decreases with iV, suggesting 
that the mean-field limit is well-defined for infinite popu- 
lations. Conversely, fixing N and increasing t, 7 appears 
to reach finite values, hence Var^ (and hence Vk) seems 
to be non-self averaging in time. This indicates that a 
mean-field description of the population dynamics might 
be appropriate for longer time-scales at increasing popu- 
lation size. 

Note that these questions are superfluous for any em- 
pirical purpose of the model, since the infinite N mean- 
field limit does not describe correctly the finite N dynam- 
ics already for a = f [T] , and the time scale of any exper- 
iment is orders of magnitude smaller than those consid- 
ered here. "Real" evolutionary time scales are longer, but 
presumably subject to changing conditions, and therefore 
beyond the scopes of this model. However, the infinite 
N limit is instructive on purely theoretical grounds, for 
understanding the main mechanisms driving the model. 

To sum up, the mean-field regime predicts that the 
infinite-population/infinitc-time speed for the mutation 
classes is Ui,- At finite N, the long-time averages for 
the number of mutations k and the associated speed Vk 
agree with the predictions of the mean-field regime, but 
are characterized by strong fluctuations. Finally, an ex- 
tended expression accounting for the finite width Varj, of 
the class f{k, t) relates Vk to the width itself, and is ver- 
ified by simulations. In this long-time regime, a number 
of assumptions of the model break down, and more com- 
plex models need to be employed. For practical purposes, 
neither the infinite population size nor the infinite time 
limit are relevant, as even in long-term experiments the 
number of generations of interest is typically sufficiently 
low. This will be discussed in a later section, when the 
model parameters are matched with empirical data from 
two laboratory evolution experiments. The next sections 
discuss in greater detail the behavior of the model at fi- 
nite N. 



C. At finite population size, the advantage and 
fitness-class distributions remain Gaussian, but their 
width is unsteady, with opposed trends. 

The large-iV long-time phenomenology described so 
far and captured in terms of mean-field estimates is not 
very relevant empirically, given the experimental range 
of times and population sizes. In the following we offer 
a more detailed summary of the behavior of finite-size 
populations at intermediate times (i.e. on the relevant 
experimental time scale « 10^ — 10^ generations, with 
{k) « 10- 10^). 

As already illustrated, the frequency distributions of 
the mutation and advantage classes are approximatively 
Gaussian for a constant advantage model. Direct simu- 
lations of the diminishing return model show that this 
still holds (Fig. |4|\) . Given the discrete nature of these 



distributions, their widths are well represented by the 
distances Lk and Lg of the foremost bin from the aver- 
age (shown in Fig.[T]A). While in a fixed advantage model 
these distances, in terms of both mutation and advantage 
classes, between the edge and the mean are constant, for 
a diminishing return model Lk is an increasing function 
of the mean number of mutation classes (fc), while Lg 
decreases with [k) (Fig. |4j3). 

The model without epistasis (and the mean-field es- 
timates) make us expect that the widths of these his- 
tograms are related to the speed of adaptation and con- 
sequently of mutation accumulation. The instantaneous 
velocities of the two distributions Vk and Vs both decrease 
with (fc). 

Qualitatively, one expects that, with changing mean 
number of mutations, the two speeds are connected by 
the decay of the advantage between k and k + \. How- 
ever, the relation between them is also affected by the 
change in the width of the distributions with increasing 
{k) . Indeed Vs decreases more rapidly and, in the long- 
time limit, even if Vk > 0, it vanishes. This implies the 
existence of a large number of sub-populations with dif- 
ferent k but essentially the same fitness, which leads to 
the effectively neutral behavior discussed in the previous 
section. Since the advantage decreases with time, and 
the distributions change width, it is interesting to ask 
under which conditions the biologically relevant regimes 
are preserved. Specifically, the condition of multiple con- 
current mutations could be lost over time. 

As discussed in the Introduction, the clonal inter- 
ference condition can be roughly defined (neglecting 
irrelevant numerical factors) by the scaling relation 
NUbln{NAsk) > 1, where Ask = sq (^(fcmax) - gi{k))) 
is the advantage of the fittest class. This condition is sat- 
isfied when the population size and the edge advantage 
Asfc are large enough, i.e. when NAs ^ 1 and NUt ^ 1. 
Considering a fixed advantage model, for typical values 
of experimental population size {N « 10^ — 10^°) and 
advantage (sq ~ 10^^), these inequalities are always sat- 
isfied, for Ub estimated in the range 10~^ — 10~^° mu- 
tations per genome per generation. However, when the 
contribution of a newly arising mutation becomes too 
small, and the advantage in fitness of the edge (Asfc) 
is smaller than 1/iV, the clonal interference hypothesis 
breaks down. Note that this implies also an establish- 
ment size larger than N. These considerations confirm 
that, while the long-time limit of this model is biologi- 
cally implausible, the experimentally relevant parameter 
regime should be captured correctly by our theoretical 
framework. 

We will now attempt to generalize the analytical es- 
timates for the finite-A^ adaptation speed available for 
the multiple mutation model to the case of diminish- 
ing returns. These "stochastic edge" estimates are based 
on the hypothesis that the only class subjected to sub- 
stantial stochastic effects is the fittest one, i.e. that 
As > Ub [ini HZ]- When As w Ub the single stochas- 
tic class approximation fails, and the dynamics is more 
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complex. Supposing that sq = 5 ■ 10~^, a = 0.02 (this 
is a much stronger epistatic effect than the one we es- 
timate from experimental data, see the following) and 
that Lk ~ 50 (of the order of the values obtained from 
simulations with this parameter set and population size 
N = 10^ — 10^"^), the advantage becomes close to a benefi- 
cial mutation rate of Ui, w 10""^ (which can be considered 
very large t22i|23j), for (fc) w 5 • 10^. This exceeds the in- 
teresting experimental range of (fc) (10^ — 10^), suggesting 
that for the relevant range we can always suppose that 
the stochastic edge approximation is valid. 



D. Analytical estimate for the adaptation speed at 
finite A^. 



Since we have verified that the advantage and muta- 
tion class histograms are both nearly Gaussian (but not 
stationary in width), it is possible to attempt a general- 
ization of the estimates applied for the standard multiple 
mutation model, Eqs. (|3| and ([5]). Supposing a slow in- 
crease of the width of the distribution with k and assum- 
ing that the width of the histogram is stable while the 
new class is being established, the mean of the distribu- 
tion moves from {k)t to {k)t + 1 during the establishment 
time Tk and Lk « Lk+i- For (k) ^ Lk > I (and (fc) not 
too high due to the condition of sufhciently large As) 
the advantage of the edge respect to the mean class is 
Ask = so(^" - (fc)") ~ soaLkk''-\ 

Thus, imposing the condition that one new mutation 
class is established at the edge of the histogram gives 



dt 



Ub 



2soaLfcfc""^ 

(2sQaLk+i{k 



1) 



(12) 



As in Eq. ([3|, the last term in this integral is the estab- 
lishment probability of a new fittest class, while the first 
is the rate of beneficial mutations from the previous fit- 
ness class, which is born with size 75 r"*" , ^--i and grows 

exponentially. Note that imposing a — 1, Eq. (12) be- 



comes Eq. ([3|), recovering the standard multiple-mutation 
model. 

An estimate of Tk can be obtained integrating the 
above expression 



1 = 



{k + lf 



Soa(Lfe-l) fc2(a-l) 



oSo(-Lfc-l)fe° Tfc 



1], (13) 



where we used the approximation Lk ~ ifc+i- Under the 
assumption that k is not too large (i.e. the advantage 
of the fittest class is sufficiently high) the contribution 
of the integration boundary t — can be neglected and 
the expression of Tk is the equivalent of Eq. Q for the 
non-epistatic case. 



Tk 



soa{Lk - l)k° 



log 



soa{Lk - l)P("-i) 
{/fc(fc + l)"-i 



(14) 



The above expression can be further simplified assum- 
ing that both k and Lk are large enough so that fc -|- 1 ~ A: 
and Lk — 1 — Lk and neglecting the logarithmic term in 
Lk leading to the following expression 



1 



Tk 



soaLkk 



— T log 

Q— 1 ^ 



soa 



Ub 



(15) 



This indicates that for intermediate k, the estimate of 
the standard multiple mutation model is valid provided 
the diminishing return advantage function sq/c" is substi- 
tuted to the constant advantage. For sufficiently large fc, 
expansion of the exponential in Eq (13) gives Tk ~ ^/Ub, 



compatible with the mean- field result. 

The time Tk is related to the instantaneous speed of 
the mutation class histogram Vk- The speed Vs, can be 
obtained knowing that during the time the mutation class 
histogram travels by one class, the advantage histogram 
has to move by the relative fitness between the newly 
added fittest class and the previous one, SQak°'~^. Using 



the simplified Eq. ( 15 ) (which neglects logarithmic terms 
in Lk), one obtains the speed 



and 



Vk 



Soak" 

Tk 



SpaLkk^"'' 
log(^ 



sla^Lkk^^''- 



log(- 



(16) 



(17) 



The second part of the estimate involves the normal- 
ization condition ([5]). Assuming (as in the standard es- 
timate) that the largest term of the histogram domi- 
nates, we need to evaluate the time necessary for the 
fittest class to become the class with mean advantage, 
whose size is order A^/2. If the fittest class has k mu- 

Its growth 



tations, its establishment size is 75 t—tt. 

will be roughly exponential, with a rate that decreases 
while it gets closer to the mean. We estimate its growth 
by its mean growth rate during the time t^. Immedi- 
ately after establishment, its relative growth rate will be 
SQCkifcfc""^, while its rate will tend to zero when it gets 
close to the mean. Thus, on average, we can assume that 
it grows exponentially with rate . 



This argument leads to the equation 



A^/2 



1 



2soaLkk°' 



(18) 



which corresponds to Eq. ([5|, and implies the equivalent 
of Eq. 



k" 



SQahk 



(19) 



In order to estimate Wj, we need to determine how 
much the histogram of fitness advantage has progressed 
during the time from the establishment of the k-ih. 
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FIG. 4: (Color online) The histograms of fitness advantage and mutation classes have nearly Gaussian forms. 
While adaptation slows down, the latter histogram expands while the former becomes increasingly peaked. 

(A) Histograms of the mutation classes (top) and advantage classes (bottom) obtained from simulations averaging over 200 
realizations at different generations (different symbols, see legend). The parabolic form in the semi-log plot indicates that they 
are approximately Gaussian. The establishment size is represented as a dashed line in the top panel. Solid lines connecting the 
symbols are guides to the eye. (B) Simulated data for the widths of the mutation class histogram Lk (green squares, top) and 
of the fitness advantage histogram, Ls (blue circles, bottom), plotted as a function of the mean number of mutations (k). The 
continuous line represents the theoretical estimates of the width (see Eq. ( 21 1). (C) Plots of the speed of mutation accumulation 
Vk (top, green squares) and of adaptation Vs (bottom, blue circles), as a function of the mean number of mutations (k). 
Continuous lines are the corresponding theoretical estimates (see Eq.(|16[) and (17l). Note that since depends logarithmically 
on Lk, the estimates for both speeds are in satisfactory agreement with the simulated data even if Lk is approximated more 
roughly. The parameters used in the simulations are N = 10^, Ut ^ G - W~^, So = 0.1, a = 0.2, compatible with those 
estimated (see Sec III El. Averages are performed over 200 realizations. 



mutation class. We assume that during this time Lk is 
roughly constant, so that after time r^, k + Lk muta- 
tions are established, and the advantage of the edge has 
reached soaLk{k + Lk)"^^- 

This allows to estimate Vg as the advantage gained di- 
vided by the time r^, i.e. 



(gpgLfc)^ k^-^k + LkY 
2 log(iVsoaLfefc" 



(20) 



Eq. (171 and (20) together allow to determine Li^ 



which can subsequently be used to obtain the speed 
of adaptation w^., or of fixed mutations Vk, using e.g. 
Eq. ( 17 1 and ( 16 ). Assuming that k+Lk — k and neglect- 



ing the logarithmic corrections in Lk, as made to obtain 
(201, we can obtain the following closed expression for 



Lk 



Lk 



21og(iVsoafc"" 



log 



(21) 



Comparison with simulated data shows that the above 
expression for Lk, despite the rough approximations 
made, is a reasonably good estimate of the width of the 
distribution (Fig.|4^). In particular, the speed of adapta- 
tion and mutation accumulation (Fig. |4]C) are well cap- 
tured by our analytical description. We should stress 
again that these expressions can be considered valid for 
intermediate values of fc, as explained above. For large 



k, the integration boundary in Eq. ( 13 ) cannot be ne- 
glected, and it can be verified that Vk = 1/Tfe tends to 
a constant, Ub with the approximations taken, restoring 
the correct mean-field limit. 



E. Parameter matching of the model with data 
from laboratory evolution experiments. 

Having explored some of the main features of the di- 
minishing return model using theoretical arguments and 
simulations, we now proceed to compare it to experimen- 
tal data. We have made very clear that this model is a 
very crude description of any realistic situation. How- 
ever, its advantage is that it is simple and based on 
few parameters, so that estimating some of the model's 
fundamental parameters from data is a relatively simple 
task, which could be instructive [T3]. 

We developed a parameter matching procedure that 
produces an estimate of the beneficial mutation rate. 
This estimate depends on the specific advantage model 
chosen, but also allows to select models according to how 
well they resemble the data. Comparing the simulated 
dynamics with the experimental one, it is possible, for 
each specific data set, to determine an optimal functional 
form of the advantage g{k), thus narrowing down the es- 
timate of Ub- 

We analyzed fitness/mutations data from two labo- 
ratory evolution experiments using the three diiferent 
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variants of the diminishing return model described in 
Appendix [Aj The power law model, where the advan- 
tage is described by SQg{k) — Sgk"', is the main case 
presented in the former sections. The advantage func- 
tions of other two models considered have a logarithmic 
{g{k) — ln(/c -I- 1)) and an exponential {g{k) ~ ^i-^q ) 
dependence on k. These functions can be derived, with 
some approximations, as partial sum of the harmonic and 
geometric series, which allows to produce simple expres- 
sions of the antagonistic effect of each added mutation. 
These two model variants share most of the qualitative 
features of the main formulation. 



1. Experiments analyzed 

We applied our model to two experiments of con- 
trolled evolution. The first are the first 20000 gener- 
ations of the well-known "iJ. coli long-term evolution" 
experiment [21 IMl EO] , while the second is a chemostat 
experiment performed by some of the authors |31j . 

The two experiments concern two bacteria (respec- 
tively Escherichia coli and Acinetobacter baylyi), were 
performed using distinct propagation techniques (serial 
dilution in batch and chemostat respectively) and also 
their duration is quite different both in terms of gener- 
ations (2 • 10^ and 3 • 10^) and experimental time (« 10 
years and ~ 4 months). Another remarkable differ- 
ence between the two experiments is the population size, 
which varies every day by two orders of magnitude (be- 
tween « 5 • 10^ — 10*) for the serial dilution experiment 
while it is very large and approximately fixed in the A. 
baylyi experiment (« 3-10^''). The two experiments share 
interesting features suggesting that a diminishing return 
model might be applicable to describe the evolutionary 
dynamics of the populations. Firstly, their duration in 
terms of generations is long enough to observe a decel- 
eration of fitness increase P^. Secondly, the large (ef- 
fective) population size suggests that the clonal interfer- 
ence regime might be relevant in these experiments. The 
simultaneous presence of different genotypes within the 
population has been verified in both experiments |311I32| . 
Moreover, the decrease of the beneficial effect of the first 
five fixed mutations in the E. coli serial dilution experi- 
ment has been recently demonstrated [T3]. The effect is 
more complex than the description of diminishing returns 
given here. However, as suggested by the authors, one 
can surmise that a simplified model including epistatic 
interactions might be useful to roughly describe this phe- 
nomenon. 

A detailed description of the data and the choices made 
for the analysis procedure is given in Appendix [B] Here 
we give a brief description of the main features of the two 
experiments. 

The A. baylyi experiment studied the population dy- 
namics in a chemostat using a minimal medium supply 
for about four months, at a dilution rate D w 0.7h^^. 
The use of chemostat allowed to grow a large population 



(TV 3-10^°) under controlled conditions for a fairly long 
time. Since the number of individuals is large, it is ex- 
pected that different sub-populations will grow in paral- 
lel in clonal interference regime. This has been confirmed 
by population sequencing data ,31j. Additionally, several 
single clones where isolated from samples collected and 
frozen at different time. 

The maximum growth rate (fJ-max) of 21 isolated clones 
and of the the original strain introduced in the chemo- 
stat (wild type) has been measured in batch, fitting the 
growth curve during the exponential phase. We consid- 
ered the maximum growth rate measured in batch as in- 
dicative of the fitness of the population into the chemo- 
stat, defining the normalized fitness of the clone i — th as 
Wexpii) — e'^'""^'*-'"'^"'^. A more detailed description of 
this procedure is given in App.jB] Note that, since the ref- 
erence fitness value is given by the ancestral growth rate, 
?i'e2:p(WT) = 1. A whole-genome sequencing was per- 
formed on two single clones isolated at the end of the ex- 
periment (AB2800b and AB2800a) and on the ancestral 
strain that served as a reference for identifying mutations 
in evolved clones. A total of 11 mutations were detected 
into the evolved clones, eight of them in common between 
the two. Additional sequencing has been performed in 
the remaining clones selected at different generations on 
PGR fragments encompassing the mutated loci identi- 
fied in the two end-point clones in order to reconstruct a 
sketch of the history (or the genealogy) of the mutation 
appearances. 

The E. coli long-term evolution experiment concerns 
twelve E. coli populations evolved in parallel for about 
5-10^ generations in batch. Serial dilution was performed 
daily, allowing sa 6.64 generations each day and an effec- 
tive population size of « 2 • 10^ individuals. We used 
the mutations and fitness data of the population desig- 
nated Ara-1 referred to the first 20000 generations, as 
given in ref. |3j. The relative log-fitness of the ances- 
tor and evolved populations(that is essentially the ratio 
between the Malthusian parameters, ipi = iXi/no) was 
measured through competition experiments (see ref. |33j ) 
every 1000 generations. We define in this case the mean 
fitness of the i~th sub-population as Wexp{i) — e^''''^^^^, 
where /i = ln(2) is the approximate mean growth rate 
(see ref. [3^ and Appendix [b|) . The fitness of the refer- 
ence strain is again Wexp = 1- Genome sequencing was 
performed on samples from generation 2000, 5000, 10000 
15000 and 20000 as well as on the ancestor. A total of 45 
mutations were found in the most evolved strain, most 
of which were stable in later clones [3]. 

For all the analyzed clones in both experiments, it is 
possible to associate fitness values with numbers of muta- 
tions, and thus bridge with the parameter k in our model. 
For the A. baylyi experiment the estimate of the number 
of accumulated mutations are inferred from the geneal- 
ogy, and thus have to be considered as a lower bound for 
all the clones, except for the ancestor and the two fully 
sequenced clones where they are measured directly. How- 
ever, the indications of k from the population sequencing 
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data are compatible with the inferred values of k 



2. Matching procedure and estimate of Ut 



Estimate of 
Fitness advantage 
s„g(k) from data 



Simulations with 
fixed g(l<) and 
varying Ub 




Time steps 



experimental <k> 



Estimate Ub 
using 
experimental time 



Experimental 
number of 
generations 



Beneficial Mutation Rate Ub 

FIG. 5: (Color online) Sketch of the parameter match- 
ing procedure leading to an estimate of the beneficial 
mutation rate from experimental fitness data. Top: 
the experimental data for fitness as a function of the number 
of mutation (blue symbols) allow to estimate the advantage 
fitness function sog{k) from a fit (continuous green line). Mid- 
dle: simulations are run using the estimate advantage func- 
tion for different values of Ub, which remains undetermined. 
This leads to different predicted dynamics for the number of 
acquired mutations (dotted, dashed, and dashed-dotted lines) 
and the fitness increase in time. These predictions can be 
matched with the experiment to estimate Ub- In particular, 
the predicted number of time steps necessary to reach the 
final experimental number of mutations varies with the ben- 
eficial mutation rate (filled circles). Bottom: The estimate 
of Ub is obtained matching the predicted final time with the 
experimental one. The best value of Ub (blue star) is obtained 
when the predicted number of time steps necessary to reach 
the final number of mutations (determined in the experiment) 
corresponds to the number of experimental generations (hori- 
zontal continuous red line). We applied this procedure to the 
three model variants for the diminishing return described in 
the main text. 

Figure [5] summarizes the parameter matching proce- 
dure. The first step finds the best-fitting parameters 
for each functional form of the fitness advantage func- 
tion g{k). This uses data on fitness values and num- 
ber of acquired mutations, using the definitions given in 
the previous section. Since the number of experimental 
points is low (5 for the E. coli and 21 for the A. baylyi 
experiment- corresponding to 9 different values of k), it 
is possible to obtain good fits with different functional 
forms of the advantage sog{k). The results of the fit us- 
ing the three different functional forms considered here 
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FIG. 6: (Color online) Estimate of the beneficial muta- 
tion rate from different experiments. Estimate of the 
beneficial mutation rate using data from Jezequel et al. [31j 
(left panels) and Barrick et al. |3j (right panels). Top: es- 
timate of the fitness advantage function sog{k) from data of 
advantage as a function of mutation number (red symbols), 
as described in the top panel of Fig. [5] The different lines 
correspond to the three model variants for the diminishing- 
return advantage (power-law dot-dashed purple lines, loga- 
rithmic long-dashed green lines, geometric blue dashed lines, 
as in legend). For sake of simplicity for experiment A. bay- 
lyi (panel Al) for each k only the mean value of the advan- 
tage is shown. Complete data are reported in ref. 31 . Mid- 
dle: estimates of Ub, from different model variants, obtained 
matching the time for which {k){texp) = {k)exp predicted by 
simulations with a given sog{k) to the experimental number 
of generations, as described in the bottom panel of Fig. p] (the 
different line styles refer to the model variants as above)! For 
A. baylyi experiment we considered a total of 2850 genera- 
tions. Bottom: comparison of the performance of the three 
model variants with the estimated values of Ub. For the Jeze- 
quel et al. data (left panel), the power law and geometric 
variants are quite close, but power law model for the dimin- 
ishing return gives the best agreement, especially consider- 
ing the data relative to the number of acquired mutations 
as a function of time (not shown). This choice leads to the 
estimate Ub ~ 10~^ — 10~^. For the Barrick et al. data, 
the logarithm and power-law models perform better, and give 
equivalent estimates of Ub ~ 10""'. 



(power law, logarithm and geometric) is shown in the top 
panel of Fig. [6j In a second step of the procedure, sim- 
ulations are repeated for a wide range of values of Ub. 
An estimate of this parameter can be obtained using the 
mean number of mutations at the end of the two exper- 
iments (Fig. [5|. In our case, (k) = 11 and (fc) = 45 for 
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the A. baylyi experiment and E. coli experiment respec- 
tively. In the model, for a fixed interval of time steps, the 
number of accumulated mutations decreases significantly 
and monotonically with Ub, and the number of time steps 
necessary to reach (fc) = 11 and (fc) — 45 depends on Ub- 
Thus, for each model there is a single value of the benefi- 
cial mutation rate that verifies (fc)(iej.p) = (fc)exp (Fig-[5] 
and middle panel of Fig. |6]). Note that we are referring 
to the estimated values of k as mean values assuming 
that the sequenced clones are representative of the pop- 
ulations. For the A. baylyi experiment, the uncertainty 
on the number of mutations present in the clones is a 
relevant source of error. On the other side, even if for 
the E. coli experiment the number of mutations is exact, 
and referred to single clones, the fitness data we associ- 
ated are population mean fitness data, that can introduce 
an error in the estimate. 

The procedure applied so far gives some estimated 
values of the beneficial mutation rate, which how- 
ever, depend on the advantage model, varying at this 
stage up to three orders of magnitude. Nevertheless, 
all the estimates are roughly in the expected biological 
range (10^ - 10^ see Ref. [221I23]). The third step of the 
estimate procedure allows to select between advantage 
model, and is based on the comparison of simulated and 
experimentally measured dynamics for the fitness and the 
number of mutations as a function of time. Note that this 
is not trivially equivalent to the fitness as a function of 
/c, but contains the effects of the population dynamics 
in presence of clonal interference, for which the model 
provides a description. Since different functional forms 
of the advantage imply a different behavior of the pop- 
ulation, one can check which of the advantage functions 
gives the closest description to the experiments. Thus, 
the experimental data of the advantage as a function of 
the time can be used to discern between power law, log- 
arithmic and geometric model, ultimately selecting an 
estimate of Ub- 

A qualitative comparison between data and simula- 
tions shows that for the A- baylyi experiment the power 
law and geometric advantage model (with parameters 
So = 0.420, a = 0.194 and sq = 0239, q = 0.631 re- 
spectively) best describe the increase of the fitness with 
time (bottom panel of Fig. [6| . Comparing the increase of 
the number of mutations as a function of time of the two 
models suggests that the power law model better resem- 
bles the data. For this reason, we speculate that, despite 
of the large errors in the fits, the power law model could 
perform better in this parameter-estimation procedure 
and the preferred estimated value of beneficial mutation 
rate is around the value C/f, w 3 • 10^^. In the E- coli 
experiment, the power-law and logarithm model result 
describe the data best, and are roughly equivalent. The 
estimate parameters are sq = 0.220, a — 0.258 for the 
power law and sq = 0.158 for the logarithm model. The 
estimates for the beneficial mutation rate are also essen- 
tially equivalent (considering the errors connected to the 
procedure) for the two preferred models [/& « 1 • 10~^ for 



the power law and Ub ^ 6 ■ 10~^ for the logarithmic ad- 
vantage model. Note that for this step of the procedure 
the data used are much more abundant, since we used all 
the fitness data measured every 1000 generations. A log- 
arithmic increase of the fitness advantage for E- coli long 
term evolution experiment is also suggested by a parallel 
work [M] . 



IV. DISCUSSION AND CONCLUSIONS 

Different laboratory evolution experiments show a de- 
crease of the fitness advantage due to newly acquired mu- 
tations and a decrease of the speed of evolution pi lTClISS] . 
This effect can be explained in different ways, and ac- 
cordingly different models have been formulated in this 
context. For example, the speed of evolution could de- 
crease because beneficial mutations with larger advan- 
tage fix sooner in the population [18] and because the 
mutation rate or the number of possible beneficial muta- 
tions decreases with time [271 HH] • Another explanation, 
suggested by recent experimental observations [HI [15], 
could be that increasing the number of accumulated mu- 
tations, epistatic interactions cause a decrease of their 
effect on the fitness. These different explanations are not 
necessary mutually exclusive, and could be stratified in 
actual laboratory evolution experiments [1]. 

We considered a simplified model, using a minimal 
number of parameters, which is a direct generalization 
of the multiple mutations model at constant advantage, 
but describes diminishing returns. Specifically, it is as- 
sumed that the selective advantage of all individuals hav- 
ing k beneficial mutations is identical, but decreases with 
k. While more complex and realistic descriptions appli- 
cable to laboratory evolution exist an advantage of the 
approach taken here is that, as the multiple mutation 
model, it depends on few parameters, and as we have 
shown, in line of principle allows matching of these pa- 
rameters with data. On theoretical grounds, this class of 
models is interesting because its dynamics is driven by 
rare events [TOl [H] , and because its behavior at finite N 
is qualitatively distinct from its mean-field behavior IJ. 

We have shown that basic phenomenology of the model 
entails a sublinear decrease of the mean number of fixed 
mutations and a steeper sublinear decrease of the mean 
advantage. This is in qualitative agreement with previous 
results using a similar model applicable in a regime where 
concurrent mutations do not occur [IB] . 

The two evolutionary speeds are related to the width of 
the distribution of coexisting advantage classes, and thus 
of coexisting mutation classes. We showed how a theoret- 
ical mean-field argument produces a relation between the 
speed of fixed mutations and the second moment of the 
histogram of mutation classes, confirmed by simulations. 
In the limit case of constant advantage, this relation re- 
duces to the one previously derived in ref. [1], related to 
Fisher's fundamental theorem. Interestingly, simulations 
indicate that for any finite N, different model realizations 
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behave increasingly differently with time in terms of both 
Vk and width of the mutation class histogram. This non- 
self-averaging property implies that even at intermediate 
times, the behavior of a realization can be quite different 
from the average. 

Within a mean-field framework, we also derived ar- 
guments showing that the mutation speed d{k) /dt ap- 
proaches the beneficial mutation rate Ub for very large 
times, which corresponds to the mean behavior observed 
in simulations at large times. In fact, as the original 
multiple mutations model, the present model does not 
include deleterious mutations. However, in the case of 
diminishing returns the assumption is more delicate, be- 
cause the advantage of fixed mutations will decrease in- 
definitely, until no beneficial mutation is able to fix. Em- 
pirically, before this happens, beneficial mutations will 
be balanced by deleterious ones [26] . This makes the (ef- 
fectively neutral) long-time limit of the model unphysi- 
cal. We have discussed how this limitation should not 
affect time-scales and parameter values that are relevant 
for current laboratory evolution experiments focused on 
adaptation [5]. A generalization of the model that in- 
cludes deleterious mutations is possible and could be the 
subject of future investigations, since it allows to for- 
mulate evolutionary questions related to the balance of 
mutations of different kinds [Ml [37] . 

Finally, for finite population size N , we were able to de- 
fine through analytical arguments the regime where the 
stochastic edge estimate of the adaptation and mutation 
speeds can be extended to the case of diminishing returns. 
In the simplest case it is possible to obtain closed expres- 
sions for the mutation class and adaptation speed and for 
the supports and of the mutation class and advan- 
tage histogram respectively. These expressions are direct 
generalizations of those obtained in the case of constant 
advantage |101 111) , and compare well with numerical sim- 
ulations. Once again, in presence of diminishing returns 
these expressions, besides the usual limitations in terms 
of parameter values, have a limited applicability in time, 
but are typically valid in the experimental range. 

The constant-advantage multiple mutation model has 
previously been applied to short-term laboratory evolu- 
tion experiments [13| . In those early stages, adaptation 
does not slow down, and the assumption of constant ad- 
vantage is justified, allowing to use the model to estimate 
empirical parameters. However, on longer but experi- 
mentally observable time scales, the assumption breaks 
down, and arguably the closest extension of that model 
that can be used to estimate parameters is the one de- 
scribed here. 

In order to show a proof-of-principle parameter- 
matching procedure for the diminishing return model, 
we considered two different experimental data sets from 
long-term evolution experiments, and defined a proce- 
dure that can be used to infer the order of magnitude 
of the beneficial mutation rate from genomic and fitness 
data, assuming the model. This procedure allows to es- 
timate the beneficial mutation rate Ub and also compare 



different functional forms of the decreasing advantage 
functions against data, through simulations of the model. 
The values obtained for the beneficial mutation rate fall 
within the range of the available measurements [551 123] j 
on the order of 10^^/10^^ mutations per genome per gen- 
eration for the E. coli long-term evolution experiment, 
and between 10^^ and 10~^ in the case of A. baylyi. 
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Appendix A: Definition of the three model variants 

This Appendix describes in additional detail the three 
model variants used in the main text. 

All variants start from the main assumption that the 
selection coefficient is dependent on the number of mu- 
tations s = Sog'{k), where g'{k) is a decreasing function 
of k, and corresponds to different specifications of g'{k). 
Every individual with k beneficial mutation has fitness 

Wk = e^^'=o = e^''^'(''\ (Al) 

and we can always arbitrarily define w{0) = 1 as the 
"wild- type" fitness. The main example of g'{k) we con- 
sidered is given by the choice of a fitness gain that de- 
pends on the number of mutations occurred fc as a power 
law. In this case, g'{k) — ak"^^ with a < 1, where the 
epistasis grows in strength decreasing a from the non- 
epistatic case of a = 1. In this case, the fitness is 

Wk = e^''=o'^°"'='°"' , (A2) 

with a < 1 for diminishing return, and a = 1 — > for no 
epistasis. 

In the general case a ^ 0, one can write 

Wk = e"""""-^-- « f,'>a{k'^'^[^+\+o(i/k)\+ai-a)) ^ (A3) 

where C, indicates the Riemann zeta function. Thus, the 
relative fitness can be expressed as 

Xk ~ e'*°('=°-<'=°» (A4) 

For the particular case a = 0, 

Wk = e^'^' = i W g^o(in(fc)+7) ^ (A5) 

which is obtained truncating the harmonic number ex- 
pansion, and where 7 = 0.57721 is the Euler-Mascheroni 
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constant. This approximation neglects the terms 0{l/k). 
Under this assumption, the relative fitness can be ex- 
pressed as 



the power law model. In the same quasi-Gaussian ap- 
proximation described for the power law and logarithmic 
model, the relative fitness is given by 



Xk 



(A6) 



These expressions show that it is essentially equivalent 
to directly assume fitness functions of the form 



Wk 



if a=^0 
, 



(A7) 



that can be expressed, neglecting terms of order 1/fc, as 
sum of powers of the number of accumulated mutations. 
In the case a = the factor 7 results a multiplicative 
factor that can be canceled since it does not affect the 
relative fitness, that is the relevant quantity for the dy- 
namics. On the other hand, this notation requires an 
appropriate rescaling of squ — > a in the case of a 7^ 0. 
Under these assumptions, the relative fitness, which ul- 
timately defines the dynamics, is identical to the one ob- 
tained considering the explicit sum of the contributions 
of single mutations, neglecting corr ections for small k. 
Moreover, the fitness functions A7 using ln{k + 1) in- 
stead of ln{k), allow to automatically include the case 
k = with the correct normalization condition w{0) = 1. 
We refer to the case a ^ as the power law model, while 
the case a — corresponds to the logarithmic model. 

In the case of an approximately Gaussian distribution 
of mutation classes with distance Lk between the fittest 
class and the class with mean log-fitness, the relative fit- 
ness of the fittest class can be expressed as 



Xk 



(A8) 



where the expansion is performed for Lk/k <ti 1 which 
should be the case in our system. In the main text, 
we use this approximation to perform the estimate of 
the adaptation speed. In the non-epistatic case {a = 1) 
the expression gives the relative advantage expression of 
the fittest class used in the standard multiple mutation 
model. Similarly, in the limit a — 



Xk 



^soLk/k 



(A9) 



The third model we use is characterized by a "geomet- 
ric dependence" of the fitness advantage on the number 
of acquired mutations. In this case the advantage func- 



tion is g'{k) 
by 



q'^ ^ with q < 1 and the fitness is given 



Wk 



(AlO) 



In other words, the advantage accumulates following a ge- 
ometric sum. As for the former models, considering only 

the final form of the advantage function g{k) = — ^ 
allows to directly satisfy the condition g{0) = 
constant factor I — q could be adsorbed in sq, as for 



l-q ' 

The 



Xfc « e " i-<J , 
and, in the case of the fittest class it becomes 



Xk 



^(k) l-q^ 



(All) 



(A12) 



Note that, while in the previous models the fitness 
is not bounded from above, the fitness for the geomet- 
ric model is finite for infinite (fc) (i.e. w{k — >■ 00) — 
1 
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FIG. A.l: (Color online)The model is effectively invari- 
ant for rescaling of the parameters. The figure shows 
the mean number of mutations ((fc), panel A) and the fit- 
ness advantage (panel B) obtained from simulations run using 
three different rescaling factors (r=l,10,100 different symbols 
as in legend). The simulated dynamics is almost unaffected 
by the rescaling procedure. Simulations are made using the 
parameters N = 10^, sq = 0.1, a = 0.2, Ub = 2 ■ 10"''. The 
data are averaged over 100 iterations, and error bars (stan- 
dard error) are smaller than symbols. 
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FIG. A. 2: The fixation probability is proportional to 
the fitness advantage. The figure shows the fixation prob- 
ability n{s) of a single clone that grows in a uniform back- 
ground having advantage s (symbols), measured from our 
simulations. The fixation probability is obtained as the ra- 
tio between the number of realizations where the beneficial 
mutator fixes and the total number of realizations (10*). The 
results give 7r(s) = 2s (continuous red line), in accordance 
with [5T]. Simulations are performed using A'^ — 10^. 
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FIG. A. 3: The relative variance of the distribution de- 
creases in time, while the relative variance of the ve- 
locity follows an increasing trend. The figure shows 
the ratio between the variance Var^ and the mean number of 
mutations (k) (left panel) and the ratio between Var„j. (cal- 
culated over different realizations) and the velocity Vk (right 
panel) as a function of time. Simulations results confirm that, 
despite Var^ increases in time, the mean number of mutations 
increases more rapidly. Since Dt < Lk, simulations results 
confirm the hypothesis Dk/{k) <C 1 used in the main text 
(see also Fig. |4|3) . The increase in variance of the velocity re- 
flects both the increase of the variance of the distribution and 
the decrease of Vk- However the error is sufficiently small not 
to affect the results on experimentally relevant time scales. 
Simulations are performed using the parameters A'^ — 10^, 



so = 0.5, a = 0.02, Ub = l- 
realizations of the process. 
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Appendix B: Experimental data and definition of 
fitness 

This section briefly summarizes the main features of 
the evolutionary experiments with A. baylyi and E. coli 
considered here, in order to describe the assumptions 
made to link genomic and growth rate measurements to 
fitness in the model. More detailed information about 
the long-term evolution experiment with E. coli, and in 
particular about fitness measurements can be found in 
refs. 1301 • Fitness data used during the parameter- 
matching procedure are reported in ref. [31 133] . A com- 
prehensive description of the experimental methods used 
for the experiments with A. baylyi can be found in 
ref. |31j . Note that in our modeling framework, as in 
most standard evolutionary models, the population size 
is kept constant and all individuals are substituted by 
newly generated offsprings at every generation, assuming 
a constant time interval between generations, although 
none of these assumptions are completely verified in the 
two experiments we considered. 



1. Acinetobacter baylyi evolution experiment 

In the A. baylyi experiment the population size is 
almost constant and, even if generations are not syn- 
chronous, the mean growth rate is held fixed by the 
dilution rate in the chemostat. Maximum growth-rate 
measurements were performed in batch on single-clone 



colonies. These values differ from the effective growth 
rate that each population can reach in chemostat be- 
cause of competition between different sub-populations 
for a limited amount of nutrients. However, the mea- 
sured values of iimax are supposed to be indicative of the 
effective fitness of the different sub-populations inside the 
chemostat. 

When a monoclonal population grows in the chemo- 
stat, its growth rate equals the dilution rate D, while in 
presence of different sub-populations the dynamics be- 
comes more complex |36j . For the sake of simplicity, and 
since the detailed dynamics is not experimentally acces- 
sible, we assume that the growth rate of all clones is 
fixed by dilution rate in the chemostat over the whole 
experiment and that the generation time is defined as 
tgen — ln{2) / D . Thus, the different maximum growth 
rates of clones measured in batch are interpreted as dif- 
ferent survival probabilities of their offspring inside the 
chemostat. 

Defining the fitness as the mean expected number of 
(surviving) offspring per generation inside the chemostat 
we assume that it is estimated by the maximum growth 
rate of mutant, i.e. w{i) = e"^' = e^"^'"'-^^^""- , where 
the index i indicates the i — th sub-population, is the 
growth rate expressed in 1/ generation, experimental data 
of l^max,i havc the units of h~^, and tgen = ln{2)/D « 
lh~^ {D K, 0.7h~^). In simple words the fitness of an 
individual is defined as the amount of offspring assuming 
it can grow by the maximum growth rate over an average 
generation defined by the chemostat dilution (and thus 
fitness is a dimensionless quantity). Since the relevant 
quantity of the dynamics is the relative fitness, during the 
parameter-matching procedure we used the normalized 
fitness {w{i) = e'^'""="-*~^""'="'"'^) to infer the functional 
form of the advantage. 



2. Escherichia coli evolution experiment 



The long-term E. coli evolution experiment is per- 
formed in batch, and serial dilution 1 : 100 is performed 
daily. The maximum population size (« 5-10^) is fixed by 
the total amount of nutrient in the medium. The number 
of doubling for each day is ^0(72 (100) ss 6.64 as derived 
from the 100- fold daily increase f3(J]. The effective pop- 
ulation size can be estimated around 2 • 10^ individuals. 
Since the speed of evolution depends logarithmically on 
the population size, the error on the estimate of Ui, given 
by this approximation should be very low. Competition 
experiments are performed between samples at different 
times and a spontaneous mutant of the ancestor, which 
is easier to track visually and has been verified to have 
almost the same fitness [33j. Relative log-fitness of pop- 
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ulation i respect to the wild type is defined as 

1 fMi\ 

In (4^) 

\f{0)wTj 

where /(I) and /(O) indicate the frequencies at the begin- 
ning and the end of the competition experiment [31 133j . 
ip{i) can be seen as the ratio between the growth rate 



of the two populations and corresponds to log-fitness in 
the model. Note again that the experimental data are 
dimensionless. The difference between the Malthusian 
parameters can be approximatively deduced as {^{i) — 
l)/i « (pi — ij-wt) where fl is the mean Malthusian 
parameter [30] . Since we are interested in expressing 
the fitness as the mean number of offspring per gen- 
eration, we use as mean Malthusian parameter fl = 
ln(2)generations~"'^. Thus we define the normalized fit- 
ness as w{i) = e('/'W-i)i"(2). 
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